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In traditional hydrodynamic theories for ionic fluids, conservation of mass and linear momentum is not 
properly taken care of. In this paper, we develop hydrodynamic theories for viscous, ionic fluids of N ionic 
species enforcing mass and momentum conservation as well as considering the size effect of the ionic 
particles. The theories developed are quasi-incompressible in that the mass-average velocity is no longer 
divergence-free whenever there exists variability in densities of the fluid components, and the theories are 
dissipative. We present several ways to derive the transport equations for the ions, which lead to different 
rates of energy dissipation. The theories can be formulated in either number densities, volume fractions or 
mass densities of the ionic fluid components. We show that the theory with the Cahn-Hilliard transport 
equation for ionic species reduces to the classical Poisson-Nernst-Planck (PNP) model with the size effect 
for ionic fluids when the densities of the fluid components are equal and the entropy of the solvent is 
neglected. It further reduces to the PNP model when the size effect is neglected. A linear stability analysis 
of the model together with two of its limits, which are the extended PNP model (EPNP defined in the text) 
and the classical PNP model (CPNP) with the finite size effect, on a constant state and a comparison among 
the three models in ID space are presented to highlight the similarity and the departure of this model from 
the EPNP and the CPNP model. 

© 2018 Published by Elsevier B.V. 


1. Introduction 

Phase field models have been used successfully to study a 
variety of multiphasic phenomena like equilibrium shapes of vesi¬ 
cle membranes [13,14], blends of polymeric liquids [2,3,17,49-51], 
multiphase fluid flows [1,19,25,34,35,38,42,54-56,58,60], dentritic 
growth in solidification, microstructure evolution [21,28,40], grain 
growth [9], crack propagation [10], morphological pattern forma¬ 
tion in thin films and on surfaces [36,45], self-assembly dynamics 
of two-phase monolayers on an elastic substrate [37], a wide 
variety of diffusive and diffusion-less solid-state phase transitions 
[11,53], dislocation modeling in microstructure, electro-migration 
and multiscale modeling [46], Multiple phase-field methods can be 
devised to study multiphase materials [54], Recently, phase field 
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models are applied to study liquid crystal drop deformation in 
another fluid, liquid films, polymer nanocomposites, biofilms and 
cells [ 18,19,25,32,34,35,38,52,54-56,58,59,61 -64], 

Comparing to other mathematical and computational technolo¬ 
gies available for studying multi-phase materials, the phase-field 
approach exhibits a clear advantage in its simplicity in model formu¬ 
lation, ease of numerical implementation, and the ability to explore 
essential interfacial physics at the interfacial regions, etc. Comput¬ 
ing the interface without explicitly tracking the interface is the most 
attractive numerical feature of this modeling and computational tech¬ 
nology. Since the pioneering work of Cahn and Hilliard in the 50’s of 
the last century, the Cahn-Hilliard equation has been the foundation 
for various phase field models [7,8], It arises naturally as a model for 
multiphase material mixtures should the entropic and mixing energy 
of the mixture system be known. 

While modeling immiscible binary fluid mixtures using phase field 
theories, one commonly uses a labeling or a phase variable (also 
known as a volume fraction or an order parameter) <f> to distinguish 
between distinct fluid phases. For instance cf> = 1 indicates one fluid 
phase while cf> = 0 denotes the other fluid phase in an immiscible 
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binary mixture. The interfacial region is tracked by 0 < <J> < 1. For 
historical more than logical reasons, most mixing energies are calcu¬ 
lated in terms of the volume fraction instead of the mass fraction in 
the literature [20,12]. Consequently, the system free energy including 
the entropic and mixing contribution has been formulated in terms 
of the volume fraction as well [20,12], given in the form F[<f>, Vcj>, ■ ■ ■]. 
A transport equation for the volume fraction <J> along with the conser¬ 
vation equation of momentum and the continuity equation constitute 
the essential part of the governing system of hydrodynamic equations 
for the binary fluid mixture, where the volume fraction serves as an 
internal variable for the fluid mixture. 

In this formulation, the material incompressibility is often identi¬ 
fied with the continuity equation 


V -v = 0. (1.1) 


This assumption is plausible and indeed consistent with the fluid 
incompressibility (1.1) only if the two fluid components in the mix¬ 
ture are either completely separated by phase boundaries when their 
densities are not equal or possibly mixed when the densities are 
identical. Otherwise, there is a potential inconsistency with the con¬ 
servation of mass as well as conservation of linear momentum. This 
inconsistency has been identified in [38], but ignored by many practi¬ 
tioners using phase field modeling technologies for hydrodynamical 
systems. We note that this inconsistency occurs only in the mixing 
region of the two incompressible fluids, where the incompressibility 
condition (1.1) is no longer valid, indicating the mixture is no longer 
incompressible despite that each fluid component participating in 
mixing is. This type of fluids is referred to as quasi-incompressible 
in [38], A systematic fix to this problem for mixtures of incom¬ 
pressible viscous fluids was given by two of the authors in [31], 
where the divergence free condition is modified to accommodate the 
quasi-incompressibility. 

In modeling of ionic fluids, one recognizes that the size of ions 
matters in most ionic solutions, in particular in the ionic solu¬ 
tions in which life occurs, in the ocean, and of course in the very 
crowded conditions found in and near electrodes in batteries and 
electrochemical cells, in and around enzymes, ionic channels, trans¬ 
porters, and nucleic acids, both DNA and RNA [65]. Ionic solutions 
are hardly ever ideal: ionic size is almost always important. In mul¬ 
tispecies ionic fluids above a certain concentration or under certain 
length scales, the size of the ions matters so that the same incon¬ 
sistency issue in the models for ionic solutions arises again. That 
is one can not simply use the solenoidal condition in the veloc¬ 
ity field as a proxy for the material incompressibility. A theory for 
multispecies ions of incompressible fluid flows that respects the 
material’s mass conservation and momentum conservation needs to 
be developed. 

This paper aims exactly at developing such a theory for a mixture 
of ionic fluid flows of multiple ionic species, in which the ionic den¬ 
sities are unmatched and different from that of the solvent, and their 
size effects are non-negligible. We require the theory to be dissipative 
while conserving mass and momentum. One targeted application of 
this theory is in ion channel modeling [ 15,1 6,26], Ion channels provide 
enough data to distinguish between theories because measurements 
are available over a wide range of conditions [5,6]. Hundreds of chan¬ 
nel types are studied every day because of their biological and clinical 
significance [65]. Concentrations and electrical potentials are con¬ 
trolled in experiments and these provide sets of values for boundary 
conditions of mathematical models. Fitting the entire set with one 
set of structural parameters allows robust solutions of the inverse 
problem [5,6] and thus allows models to be distinguished. Other 
applications of the model include electrolyte fluids, biological fluids 
with charged bio-species, etc. This theory will be consistent with the 


mass and momentum conservation and demonstrates energy dissi¬ 
pation. In principle, a variety of transport equations can be developed 
for the ionic species should one knows the system’s energy dissipa¬ 
tion rate. In this paper, we propose two types of transport equations 
based on a generalized Onsager principle [57]. These two choices 
yield two types of species transport equations and corresponding 
energy dissipation rates. Their relations with respect to the existing 
electrolyte fluid models will be discussed in the text in details. 

The derivation follows the generalized Onsager principle approach 
[31,57], leading to two types of transport equations for each ionic 
species in the form of Cahn-Hilliad and Allen-Cahn type equations, 
respectively. Apparently, these correspond to two distinct energy 
dissipation rates. Their applicability to real material systems can only 
be confirmed if one could measure the systems’ energy dissipation 
rates. However, such measurements have not yet been made, as far as 
we know. So in most cases, people adopt one particular formulation 
over the others simply based on the leap of faith. 

For the new model, together with its limits in the extended 
Poisson-Nernst-Planck (EPNP) and the classical PNP with the size 
effect (CPNP), we will study their linearized stability on constant 
steady states. Instability of the PNP class of models is of direct biolog¬ 
ical interest. Actual biological channels invariably produce unstable 
currents [41 ] that switch ‘ instantaneously’ between open and closed 
levels in a random telegraph process called single channel gating 
[24]. Instability in the models of this paper may turn into gating 
when the models are extended to include noise sources and are 
focused on the behavior of just one channel protein. However, we 
will not pursue the complicated issue in this paper; instead, we 
will focus on introducing the modeling framework and presenting a 
set of thermodynamically and hydrodynamically consistent theories, 
and discuss their predictions in a simple 1-D case to highlight the 
departure of several previously used PNP type models from the new 
model. 

The paper is organized as follows. First we present the mathe¬ 
matical formulation of hydrodynamic phase field theories for mul¬ 
tispecies ionic fluid flows and various plausible formulations of the 
transport equations giving rise to the total energy dissipation. Then, 
we examine the theory in ID geometry to compare the theory 
with some existing PNP models with and without the size effect 
[15,16,26]. Finally, we provide a concluding remark. 


2. Quasi-incompressible hydrodynamic models for ionic fluids 

We develop hydrodynamic models for a viscous, multispecies 
ionic fluid in an isothermal condition, in which mass, momentum 
conservation and the total free energy dissipation are preserved. The 
governing system of equations in the model includes the transport 
equations for all the ions, the Poisson equation for the electric poten¬ 
tial, and the conservation equation for mass and linear momentum 
of the fluid, respectively. 

2.1. Mass and momentum conservation equations 

We first present the mass and momentum conservation equation. 
We consider the transport of viscous, ionic fluids made up of N dif¬ 
ferent ionic species, each of which consists of a type of ionic particles 
of the identical size. Here, we tacitly assume the viscous solvent par¬ 
ticle is a type of ions with a zero charge [4,29,30,47]. We denote the 
number density for each type of ions by n,-, i = 1, ■■■ ,N. The elec¬ 
tric potential generated by these ionic particles is denoted by <t>. We 
denote the volume of each individual ionic particle by v,- and the 
mass by m, for i = 1, ■ ■ ■ ,N, respectively. Then, there is a constraint 
2)jLi n i v i = 1. which states that the excluded volume of the ions is 
a constant before and after the mixing. We identify i = a as the 
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solvent component which is neutral. The total density of the mixture 
is defined by 


the charge source with homogenous boundary condition and it can 
be expressed by using the Green’s function G(x,x') as 


N 

p = '^m i n i . (2.1) 

i —1 

We denote the intrinsic density of the ith species by p,- = m,/v,' t 
which is a constant. Then, it follows that 


t hn( x ) = -J (g (x,x 0 (e 0 (x') + X, z i en i (x'))) dx'. (2.6) 

Then the variation of the electrical energy F e = f a p e (j't’n + 4> e ) 
with the ion density n,- is 


N N 

P = X Pi n,V * = X foPi’ ( 2 - 2 ) 

i —1 i=l 

where <£, = nM is the volume fraction of the ith ion in the mixture. 
We introduce the mass averaged velocity v. Then, the total mass and 
the linear momentum conservation yield 


— = z,e<(>„ + Zje4>e = Zje<t. 
brij 


The equation for the total electric potential is 

V-(eV<t>) = - (e 0 +; XiZierii ), 
d’lan = i’o(9fl)- 


(2.7) 


( 2 . 8 ) 


+ V^PV) = 0, 

p^|^ + v-W^ = V-r + F (e) , (2.3) 

where r = -p 0 1 + r v is the total stress tensor, p 0 is the hydrostatic 
pressure, r v is the extra stress tensor and is the interfacial force 
that yields the Ericksen stress for the mixture fluid system. We next 
turn to the derivation of the transport equations for the ions. 

2.2. Transport equations for the ions 

The free energy of the system is prescribed as F = / n /[n lt 
■ ■ ■ ,nw]dx, where Cl is the material volume, and the density of the 
free energy functional is defined by [43,44] 

/["l,- - .nw] = fe s r X!li ^( ,nn i- 1 ) + P e Q <I, n + 

+ J K{^-y)G(ln i }f =1 {x)Ani]f =l (.yfjdy, (2.4) 

where k B is the Boltzmann constant, T is the absolute temperature, 
Nj is a generalized polymerization index for the ith ionic particle 
(N a = 1), p e = e 0 + YaL i A en i is the total charge density, z,- is the 
valence for type i ion and also denotes its sign (for solvent, we note 
that z a = 0), e is the unit charge, eo is the permanent charge density 
in the system, 4> n is the electric potential generated by the total 
charge, <t> e is a given external electric potential which is independent 
of the total charge and the total electric potential is <J> = <t>„ + <t> e . 
The first group in the sum represents the entropic contribution of the 
ionic particles to the free energy, the second part gives the electrical 
energy density of the system, and the third part gives the interaction 
of the excluded volume effect and the long-range interaction among 
the ions of finite sizes. 

The electrical energy density in the given external electric field is 
p e <T> e and in the electric field generated by the charges is lp e <t>The 
equations for the electric potentials <t>„ and «t> e are 


V-(eVcfn) = - (e 0 + XiZieni), ^ V-(£Vcb e ) = 0, 

d'nlan = 0, ‘t’elan = < t > 0 (9ff), 

(2.5) 

where £ is the dielectric constant, <t> 0 is a given boundary function. 
Here the boundary condition is Dirichlet BC, it can be changed to 
other type boundary conditions. The external electric potential <t> e is 
determined by the boundary condition with zero charge source. If 
<t>o = 0, there is no external electric potential. <t>„ is determined by 


The third part of the free energy density can be approximated via 
expansions in a differential form 


J I((x - y)G ({riilXi (x), (n,lf =1 (y)) dy « g[m, ■ ■ ■ ,n N ] 

= g(fn i lf =1 ,{Vn j )^ =1 ). 


(2.9) 


One specific form of the function g accounting for the size effect of 
the ions is given by 


g = k B T 


X' 


I fctij 'Vj -j 

H>i + X i=1 f 'I Vn 'll 2 


( 2 . 10 ) 


where the coefficient matrix E,y is symmetry. The first part in the 
energy density represents a repulsive interaction due to the finite 
size effect while the second part is the conformation entropy asso¬ 
ciated with the heterogeneous distribution of the ions in space. This 
approximate function represents the lowest order approximation to 
the interaction potential with the long-range interaction, for which 
we will adopt in the rest of the paper. The chemical potential for the 
ith ionic particle is then given by 


SF . _ 


Ni 


(In rij) + XA»"; - 7i v2|1 i 


ez,4>. 


( 2 . 11 ) 


Assuming there is no annihilation of charges between positive 
and negative ionic particles, each species’ charge and the total charge 
in the system is supposed to be conserved under the flux free 
boundary condition, 


njdx = Q,i = 1, - ■ ■ ,N, 



const, 


( 2 . 12 ) 


where Q, i = 1.---.N and C are constants and C = 0 is called 
charge neutral. Indeed, annihilation can occur in biological systems 
and ordinary bulk ionic solutions when weak acids and bases (like 
acetic acid, i.e„ vinegar, or sodium bicarbonate, i.e„ baking soda) 
are involved as components of the solution or as side chains of 
the protein that forms the ion channel. Such effects are significant 
in some cases, but they form a separate field of investigation, in 
theory, experiment, and indeed in medical practice, where they are 
particularly important. In this paper, we ignore those effects. 

We propose the transport equation for the ith ion as follows 


3n ; 

Ik 


+ V-(vn,) = Bj,i = !,■■ ,N, 


(2.13) 
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where B, is going to be determined from the total free energy dissipa¬ 
tion in the following. We note that there are two constraints of B; as 
follows, due to the constraint of n f and the total mass conservation, 
respectively. Using XtLi n i v i = XtLi = 1, we have 


dt 


( N \ / N \ N 

X n ' v ‘J + v • [Z vn < v iJ = X B < v ‘- 


(2.14) 


It implies that 


Z n „ x-,Jv „ rrij 
BjV j = V Bj —. 

i—1 1 ' ^—‘i=\ ' Pj 


(2.15) 


This gives us the first constraint on the BJs. 

In addition, from the total mass conservation and p = m i n n 
we obtain 


X rmBi = 0. 
1 = 1 


(2.16) 


This yields the second constraint on the B's. The constraints war¬ 
rant that the transport equations for each species are not completely 
independent. We next discuss two distinct ways to derive the trans¬ 
port equations for the ions and solvent following the generalized 
Onsager principle [57]. 

2.3. Formulation 1 

We denote the ath component (the solvent component) as the 
non-vanishing component in the mixture and then it follows from 
Eq. (2.16) 


B a — — 


1 


(2.17) 


The total free energy E = || v|| 2 ) dx + F of the system consists of 

two parts: the kinetic energy and the Helmholtz free energy F. Now, 
we compute the total free energy dissipation rate as follows: 


dE _d 
dt ~ dt 


I v || 2 +f dx 


-/ 

Jsi 

-/ 

Jn, 

-/ 

Jn 


*7 Wei 9nj"| , f 9/ dtij \ , c 

Vv : r- v*F lej - > Pi-*- 1 dx + / n* ( > ^^ | dS 

at J Jdn V^' =1 avn,- dt J 


Vv : r- v-F (c) + X,-! M( v-Vn i + v-Vrij) - X i= 1 W B .jdx 

Vv : Tv + X N , [(-p) m i (— - — ) -M + — Pnls, ] dx, 
^> =1 L VPi Pa/ J ) 


(2.18) 


where 3fl is the surface of the material volume fl,n is the unit 
external normal, the elastic force is identified as follows 


F (e) = XM v n„ 

1=1 

and the total pressure is given by 


P = P0- Xw n i- 
1=1 


(2.19) 


( 2 . 20 ) 


In the last step, constraint Eq. (2.17) is used. We also set the boundary 
condition 


9/ 

3 V n,- 


= 0, 


( 2 . 21 ) 


so that the surface integration is zero, i.e„ Ya=i fan 11 ' imr = 0. 

Next, we identify two forms of B, following the generalized 
Onsager principle to warrant energy dissipation of the system [57], 
They are associated with two famous transport equations: the Cahn- 
Hilliard and the Allen-Cahn equation, respectively. 

2.3.1. Cahn-Flilliard dynamics 

In the first case, we choose B ; as follows 

B, = - XL, V - A * V [(-P) m /< (k - ) ~Pk + . for i 5^ a, 

L \Pk Pa) m a 

( 2 . 22 ) 

where the mobility coefficient matrix (\y, i,j a) is symmetric and 
nonnegative definite. Then, using integration by parts, the energy 
dissipation rate is given by 


dE 

dt 


/ 1 1 \ mi; 

- — )-Pk + —Pa 
\Pk Pa) m a 


m a 


•hik 


dx + surface term < 0 (2.23) 


provided Vv : r v > 0 and the surface term is zero. For viscous fluids, 
the viscous stress tensor is given by 


Tv = 2 ?] 


D- |tr(D)I 


■ Ptr(D)I, 


(2.24) 


where D = \ (Vv + Vv T ) is the strain rate tensor, I is the identity 
tensor, q is the shear viscosity and v is the bulk viscosity. Then Vv : 
t v = 2qD : D + (it — |q^ (tr(D)) 2 > 0 is satisfied so long as q > 0 
and v - §V> 0 . The zero surface term is warranted by the following 
no-flux boundary condition: 

■t-^n , T, , / 1 1 \ mi 


= 0. (2.25) 


We summarize the governing system of equations in this model 
in the following: 


9 Mi 

It 


+ v-(v"<)= -zL (p, - ^ • 

for i 7 ^ a, 

v - v =~z^ m <(t~k) v ^ v b )mt (k ~ kY-■ + S4 ■ 

pd ii = v ■ [- ( p+ zT=, Mni ) 1 + r »]+z!l, M vni 

= V-(-pI + r„)-XH, n i v M. ( 2 . 26 ) 


and the equation for the electric potential is 
V-(eV<t>) = - (e 0 + Xi z ' eni ) ■ 


(2.27) 


where e is the dielectric constant. This model is not incompressible 
since V -v ^ 0 when densities are not identical. It is known as the 
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quasi-incompressible model [67]. This model is different from the 
previous models for ionic fluids. 

We remark that the previous models for ionic fluids assume 
the incompressible condition V -v = 0. This is valid only when 
Pi = Pj,i,j = 1, • ,JV. In this case, we end up with a self-consistent 
model as follows: 


g"' + V-(vnj) = V-A./Wlp,, -p«], for a. 

V-v = 0, 

dv , , 

p dt = v-c-pi + ^-E^w 

V -(eVO) = - (e 0 + Zjen,) . (2.28) 


In this model, the energy dissipation rate is given by 


The boundary condition for this equation system is Eq. (2.21). The 
energy dissipation rate is given by the following 

■sr^ N r, , /1 1 \ m,- 

Vv: Tv+ z^, u - 

-rt ( + ^Pa}ldx<°, (2.34) 
\ Pk Pa / Wa J J 



x\ ik 


provided (A,j) > 0. 

in the Allen-Cahn model, the charge conservation imposed by 
Eq. (2.12) may not be upheld. In order to impose the constraint 
approximately, we have to augment the free energy by adding a 
penalizing term 

u Z (J n n < - c ') 2 + L i(f n Z z < n < dx - c ) • ( 2 - 35 ) 


dE 

dt 



VV : T v + £ V[ tt -Pal-XftV [p k ~Pa] 
i,k^a 


< o. 


(2.29) 


where Li - 2 are large positive numbers. An alternative approach is to 
enforce the constraints directly by using Lagrange multipliers in the 
free energy, 


For the above two model equation systems, the following 
boundary conditions are used: 


t(L-- c ') +k {L P mdx - C } 


(2.36) 


9/ 

9 Vn, 


= o, 




m k 

Pk + — Pa 

///a 


= 0. (2.30) 


Together, they warrant that there is no boundary contribution to 
the energy dissipation and the constraints on the charge conserva¬ 
tion in the system imposed by Eq. (2.12) are satisfied. The boundary 
condition for the electric potential is the Dirichlet boundary condi¬ 
tion which is equal to a specified surface potential, and the boundary 
condition for the velocity field is the no slip boundary condition. 


2.3.2. Allen-Cahn dynamics 

Alternatively, we choose B, as follows 


,N 

B ‘ = Z*=, Xik 




rrib 

Pk+ —Pa 
lll a 


for i a, 

(2.31) 


where A,* is the mobility coefficient, we obtain an Allen-Cahn type 
transport equation for the ith ion 

§ + V -(vn,) = Zl t X ik [(-p)mk (1 - Z) -p* + ^Pa] , for i * a. 

(2.32) 


where are two Lagrange multipliers. These are common 
practices when one uses Allen-Cahn model to study multiphase fluid 
dynamics. We note that their physical validity is not widely accepted 
in the research community though. 

Note that Allen-Cahn and Cahn-Hilliard equations represent two 
different types of transport for scalar phase variables in a dissipative 
system [39]. Higher order transport equations are also possible, but 
are rarely used. Thus, we will not pursue them in this study. 

2.4. Formulation 2 


By using constraint Eq. (2.17), we rewrite the energy dissipation 
rate as follows 


dE 

dt 


-/ 

Ja 


Vv : r v 


Vv : r v 


V 

-4-<i=l 

■ v N 

£—‘i=\ 


, ^ m i 

(-P)--Mi 

Pi 


dx 


-P)—f - Pi -Lmj 
Pi 


dx, 


(2.37) 


where L is a Lagrange multiplier, which is a function of the space and 
time. If we adopt the Cahn-Hilliard equation for the ionic species, the 
right hand term B; is chosen as 


% = - Z J=1 v -v 


-P) -fJj-Lmj 
Pi 


,i = !,••• ,N, (2.38) 


where is the mobility coefficient matrix. The constraint 
Xili ntjB, = 0 implies 


The other equations are given by 


m k , 




? ' v - - ”*(s - it) h””* (s * it) 

dv f / \ 1 ,-^N 

P dt = V ‘ [- + Z i= i M n iJ 1 + TvJ + X 1= i Pi Vn > = V • (-P 1 + T v) 

v i N 

-Z i=1 n ‘ v H- 

V.(evct) = - (e 0 + X,zieni) ■ (2.33) 


Z V-\V (-p)-l-pj-Lmj 

y=i L Pi 


mt = 0 . 


It yields an elliptic equation for the Lagrange multiplier L: 


Z N ^—i N 

mimjV-XnVL = > 
«=i ’ 1 11 


V-AyV 


P] 


(2.39) 


(2.40) 


The Lagrange multiplier L is a solution of the elliptic equation. If the 
coefficient is a positive definite matrix, L is solvable in principle. In 
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a special case where A.y are constants, the Poisson equation can be 
rewritten into 


(2.41) 


Here, we don’t need to know the specific solution form for L. Then 
we have 



N 

1 N 


v 2 l = 

X X V m i m j 

M 

< 

> 

( - p)~ - Pi 

L n J 


Jj = 1 

ij= 1 


'b'S -4 


“-zLi v ' XikGk - 

The flux terms G k are given by 

\Pk Zj=i w i / 
k = 1,2.JV. 


Sij=1 hjj m i m j 


VPk 


ZjLi w j m k V[jj/m j 


SjLi w,- 


(2.43) 


The terms Wj = Xi=i = 1,2.N act as weighting 

factors. The difference between this model and the model derived in 
formulation 1 is that the correction factors are the weighted average 
terms. 

In a dilute solution, the solvent density is much larger than the 
other components, that is n a » rij for j i=- a. If we assume the 
mobility parameters Ay ~ Ajn.-dy, where A,- is a constant, then Wj = 
XI i Kj m i m j ~ A jUjinj■ Thus w a » Wj for j i=- a when mj and m a 
are not far apart, and this formulation reduces to the Cahn-Hilliard 
model derived in the previous subsection because 


ZjLj Wjin k /pj _ m k ZjLi w jm k 'Vpj/m j _ m k ^ ^ _ n 

ZjLlWj Pa ZjLlWj 


(2.44) 


For the solvent component, the governing equation of the density 
n„ is 


dn a 

dt 


+ V -(vita) = - V-AacGa ss 0. 


(2.45) 


Then we can drop the equation of the solvent component in our 
system and instead only consider the ionic components in this 
formulation. 

If we adopt the Allen-Cahn equation, the B, is chosen as follows 


The transport equation for the ith ion is given by 


3 n; 

St 


(vn f ) = y N aJ(-p)^ - 

Z^ =1 e n = iW . ; 


-Mk 


Zjii w j m kPj/mj 

ZjLi w,- 


(2.48) 


Using the same argument, if we assume the mobility parameters 

n.m. ~ X.n.m2 Thuc ia; ns. i/i/. fr»r 


5«]] 

Ay ~ Ajtiidy, then Wj = 
j a, which implies 

Et 


Zj =t Wjirifc/pj 

2-j= 

(2.42) 

ZjLiWj P«’ 



... 

-j=i w j 


-Jkxi 


and the governing equation of the solvent density n a is 


^ + V.(vn c 




( m a m a \ 

— - — -Pa+Pa 
Pa Pa J 


= 0 . 


(2.49) 


(2.50) 


This formulation reduces to the Allen-Cahn model derived in 
formulation 1. 

If (Ay) is a dense matrix, the two formulations are apparently 
different. However, if Ay = A6y, the Cahn-Hilliard equation derived 
in formulation 2 reduces to 


3 rij 

St 


V-(vrij ) = -AV 2 


+ - 


1 


-p m ±+ P 
Pk Ef=i mj 

^-,N 

X i=1 mitrikPi 


mfm k 

S-,i 


Zf=i mj 

If m, = m,i = 1, ■■■ ,N, it further reduces to 
3 n. 


dt 


+ v-(vn,) = AV 2 


1 


Likewise, the Allen-Cahn equation reduces to 
V-(vn;) = -A 


3n,- 

~3F 


1 -r—,N 

^-jvS 1=1 rt 


(2.51) 


(2.52) 


(2.53) 


Both of these have been used by some researchers in the past to 
describe multiphase materials [38], 

Apparently, formulation 2 is different from formulation 1 and it 
seems to be a more general way of deriving the transport equations 
for the ionic species. However, if we choose L such that 


-P-- Pa-lm a = 0 

Pa 


and redefine 


(2.54) 


B,= 


=xx 


III 

- p) » 


- pj - Lrrij 


(2.46) 


Ba 


1 

m 


- x B ‘ m ‘- 


ijt=a 


(2.55) 


where Ay is the mobility coefficients. The constraint X=i m;B ; = 0 
implies 2«=i [(-P)^ - Pj ~ Lm j\ mj = 0. Thus, the Lagrange 

multiplier L can be solved as follows 


L = 


Z iJ=1 ^mpnj X ij=1 *ii 


(-P)y-Pj 

P] 


(2.47) 


we recover the model derived using formulation 1. This means that 
the transport equation for n a defined in reformulation 2 must be 
modified in order to recover the transport equation in formulation 1. 
However, this modification has no impact whatsoever on the energy 
dissipation rate. 

Another remark that we would like to make on these models is 
that each model yields an energy dissipation of its own. The choice 
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of the model should therefore be made based on which energy 
dissipation rate best fits the real system to be modeled. 


is chosen as the mass density of water here. Then, we denote the 
corresponding volume scale as v 0 = Ig, mass scale m 0 = p 3 v 0 . The 
dimensionless variables are defined as follows: 


2.5. Model reformulation and reduction to existing models for 
multispecies ionic fluids 

The above models are formulated using number densities of the 
components in the fluid mixture. We can reformulate the model 
using the volume fraction 4>i or the mass fraction c,- since they 
are functions of the number density functions, <f>j = n(Vi, c,- = 
m '"', i = l, - ,N, where v, and m ; are constants, denoting the 
volume and the mass of each individual ionic particle, respectively. 

If Pi = p 0 ,i = 1,V-v = 0 and, in addition, we remove 
the entropic contribution of the solvent to the fluid mixture, i.e., we 
drop n„(ln n a - n a ), where a corresponds to the solvent component, 
from the free energy, the model reduces to the existing PNP model 
with the finite size effect [26, 27, 33]. So, all the previous ionic fluid 
models can be regarded as the model applied to the case where all 
ions are of the same mass density and the solvent effect to the free 
energy is neglected. 

Next, we compare the new model with some of its limits and 
some existing models. 


3. Binary ionic fluid model 

We consider a mixture of two distinctive ionic components 
(N = 3), where a = 3 corresponds to the solvent component, 
known as the binary ionic fluid model. The other two components 
in the fluid mixture are cations (positive ions) and anions (negative 
ions). We adopt the Cahn-Hilliard dynamics for the transport of ions. 
The governing system of equations is given by 


3n, 

~3F 


+ V-(vrij) = -V-A.,n,V I"(-p)m ( (-j- - -M -pt + ^-p 3 l, i = 1,2, 
L \Pi P3/ m 3 J 

V,v= - m (b - ti) v ’ MiV [ ( - p)m ‘ (b ~ b) + SH • 


dv 


P -fa = V-(-pI + r v )-X n < v M. 


V-(eV4>) = - 


^e 0 + X z ' en ^’ 


(3.1) 


Here, we assume the mobility matrix is \,j = the mobility of 

each ion is only dependent on its own number density. The spatial 
gradients of the chemical potentials are given by 


v Hi = j^ Vn i + eziVcb + k B T [$ 12 Vnz + SnVrq-yjVV 2 *!,], 
Vn 2 + ez 2 V<t> + + § 22 Vn 2 - y 2 VV 2 n 2 j, 


Vp 2 = 


N\n i 
k B T 
N 2 n 2 


Vp 3 =k B T- V ^ n| - V2Vn2 


■ V2H2 - Vi H] 


(3.2) 


no 


t = 


lo' 


- to - ^ - f o 

v = r v, p = -^p, Pi = u Pi. 
‘0 Polo m 0 /o 

(3.3) 


Then, the dimensionless parameters are given by 


Pi 


mo , 


rj =—, rf = — ,i = 0,1,2, p= —, pi =—, 

1 v 3 ' m 3 Po Po t 0 1 

-to - to 

T] = —jT], V = —j V, 

Pol 0 Polo 

I<b = ~%k B , ly = n 0 itj, y t = ^7;, <t> = -^4>, e 0 = 

molg % m 0 l 2 en 0 


m 0 

Zf = zi, e = - 7 ^—.e. 

S 2 tgU 0 


(3.4) 


We set = 1 to obtain t 0 = J and also set n 0 v 0 = 1 to obtain 
n 0 = 2_. It’s easy to find that r™ = pirj for i = 1,2 and r^ 1 = rJJ. 
For simplicity, we drop the ~ on the dimensionless variables and the 
parameters. The system of governing equations for the binary ionic 
fluid model in these dimensionless variables are given by 


+ v-(vui) = -V.MWHiP-M- + Ot 3 ], i = 1,2, 

V-v = - ^ =i RjV• \,njV[ —Rjp - pi + r™p 3 ] = ^ =l Rf 

dv . . ^ i3 

P^ = V.(-pI + T,)-^ = 1 n,V W , 

V-(eVt>) = - 


3 n,- 


(eo + Zi =1 z > n ')’ 


+ V-(VUi) 


(3.5) 


where the parameters R, = ( i 1 


('■o m )( p1 ' 


for i = 1,2, the total 


mass density p = 1 — Ritit — R 2 n 2 , the solvent’s number density 
n 3 = - r!(ni - r^n 2 . The spatial gradients of the chemical potentials 

are 


Vpi = 


1 


N,n, 

1 


Vn 3 + Zi V<t> + §i 2 Vn 2 + Vtlj ~7 3 VV 2 ni, 


vp 2 = Or Vn 2 + Z 2 Vd> + § 1 2 Vni + fo 2 Vu 2 - y 2 VV 2 n 2 , 
Vp 3 = 


N 2 n 2 
-rj'Vnj - r^Vn 2 


n 3 


(3.6) 


in the following, we refer to the model as the full model, where the 
word “full" means that the model respects all conservation laws and 
accounts for the finite size effect and the solvent entropy. 


3.2. Models at regimes of two distinct length scales 


Where we assume that § 3i = = y 3 = 0, i.e., the interaction 

between the ions is dominant. The entropic contribution only shows 
up in the chemical potential of solvent (p 3 ). 

3.1. Nondimensionalization 

We use a characteristic time scale to, length scale lo, and mass 
density scale p 0 = p 3 , and the characteristic number density n 0 to 
non-dimensionalize the physical variables. The mass density scale 


We examine the dimensionless full model at two distinct length 
scales. If we choose the length scale l 0 = 10 ~ 9 m = Inm, we have 
the time scale to = 1.55 x 10~ n s. If we choose the length scale 
lo = 10 _7 m = lOOnm, we have the time scale t 0 = 1.55 x 10 _6 s. 

We set the first type ion is the positive ion with valence Z\ = +1 
and polymerization index tV 3 = 1; the second type ion is the neg¬ 
ative ion with valence z 2 = -1 and polymerization index N 2 = 1. 
The values of the ratios of the ions’ volume, mass and density are 
tabulated in Table 1. The density ratio of the solvent and two ions is 
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Table 1 

The ratios of volume, mass and density. 


Ratios 


Pi 

pi 

r v r v 

'l '2 

r m ytn 

'l '2 

Values 


0.5 

2 

2 1 

1 2 

Ratios 


r v = 
'o 

. r m 
'o 

Ri 

Ri 

Values (/ 0 = 

1 nm) 

40 


0.025 

-0.025 

Values (/ 0 = 

lOOnm) 

4 x 

10 7 

2.5 x 10" 8 

-2.5 x 10- 8 


p 3 : Pi : p 2 = 1 : 0.5 : 2, the volume ratio is v 3 : v 3 : v 2 = 1 : 2 : 1. 
The size differences of the three components are distinct. The param¬ 
eters R\,R 2 are 0(10 -2 ) in the smaller length scale lo = inm. The 
compressibility of the flow (V -v 0) in the full model cannot be 
neglected. 

If the densities of the three components are the same, i.e. p 3 : 
P! : p 2 = 1 : 1 : 1, we have R 3 = R 2 = 0, the flow becomes 
incompressible. Furthermore, when the density differences are dis¬ 
tinct, but the larger characteristic length scale l 0 = lOOnmisused 
in the dimensionless system, the values of parameters R\,R 2 are very 
small, as in Table 1. If the corresponding terms of R t are dropped 
from the full model, the model reduces to a model that we call the 
extended PNP model (EPNP), in which the flow is incompressible: 


-gp + V-(vnj) = —V-XjttjV [-w + r™ft 3 ], i = 1,2, 

V-v = 0, 

p w = v ' ( ~ pI+Tv) - xLi n,VM ’ 

V-(eV<t>) = - (e 0 + mi') . (3.7) 


Here the mass density, the stress tensor and the chemical potentials 
are the same to those in the Full model. 

Furthermore, if we neglect the entropic contribution of the 
solvent to the fluid mixture, i.e., we drop n a (lnn a - l),a = 3, from 
the free energy, we get the classical PNP model with the finite size 
effect (CPNP) [15, 16, 22, 23, 48]. This is equivalent to removing the 
p 3 terms from the equations of the above EPNP model. To be clear, we 
note that the commonly used classical PNP model does not include 
the finite size effect. 

When the characteristic length scale used is l c = Inm, the 
parameters R\,R 2 are not small. So, the full model must be used. The 
model is indeed different from the limiting PNP models even with 
the finite size effect. Note that this is the length scale regime that is 
applicable to the ion channel problem. The other model parameters 
are listed in Table 2. 


3.3. Comparison of the full model with the limiting PNP models in 3D 
space 

We compare the full model with the EPNP and CPNP models in 1D 
space, assuming the system is homogeneous in the ( y,z ) directions 
and depends only on x and time t (i.e., the variables are functions of 
(t,x).) The domain for x is assumed finite given by Cl = [0,L X [. The 
governing equations of full model in ID are given explicitly by 


-jf + (Vini)' = [-lh(p)'-(|Ui)' + ff(p 3 )']}', 

^ + ( v i n 2 )' = -{A 2 "2 [-R 2 (p)' - (M 2 )' + b?(p 3 )']}'. 

(vi)' = Ri + (vini)'j + R 2 [^ +(vin 2 )'j, 

+ v i( v i)') = -CP)' + (|h + p)( v i)"-[ n i(Pi)' + " 2 (^ 2 )' + n 3 (P3)'], 

P (lJF +Vl(V2) ') = b( V 2 )"> 

p (lJF +Vl(V3) ') “ b( v 3 )"> 


(O)" = - 


2 

e 0 + X z i"i / £ ’ 

i= 1 


(3.8) 


where (•)' = ^r, (•)" = an d the gradients of the chemical 
potentials are 


(rt)' 

(mJ 

(Ms)' 


^l-(ni)' +zi(<h)' + g 12 (n 2 y + €n(ni y-7,(n,r, 

j^(n2)' + z 2 W + & 2 (ni)' + £,22(n 2 y - y 2 (n 2 )'", 

-r\{ n,y-r v 2 (n 2 y 
n 3 


(3.9) 


The unknowns are n 3 , n 2 ,p, v i, v 2 , v 3 , <t>, which are fully coupled. 

The ID EPNP model is much simpler now. First, from the incom¬ 
pressible condition (Vi)' = 0, we find that V] = 0 due to the fixed 
boundary condition of v. Then, the pressure p and v 2 ,v 3 are deter¬ 
mined from the momentum equation. The independent unknowns 
in the EPNP model are then ni,n 2 ,<J> and the ID governing equations 
are given by 


dn, 

~df 

dn 2 

~df 

(*)" 


-{\ini [-(Mi)' + )']}'. 

-{\ 2 n 2 [-(M 2 )'+ r™(M3)']}', 


eo 




/£, 


(3.10) 


where the gradients of the chemical potentials are given by Eq. (3.9). 


Table 2 

Model parameters. 


Symbol 

Parameter 

Value (unit) 

lo = lnm 

lo = lOOnm 

V 

Shear viscosity 

1 x 10 -3 kgm _1 s -1 

15.54 

155.4 

V 

Bulk viscosity 

2.75 x 10“ 3 kgm~'s ~ 1 

42.74 

427.4 

e 

Dielectric constant 

7.08 x 10 _1 °Fm _1 

0.1145 

11.45 

cb 

Electric potential 

IV 

38.65 

38.65 

7i 

High order diffusion of 1st ion 

1.6606 x 10" 27 m 5 mol~' 

to- 4 

to- 14 

72 

High order diffusion of 2nd ion 

1.6606 x 10" 27 m 5 mol~' 

10" 4 

10- 14 

Ri 

Mobility of 1st ion 

3.1083 x 10 11 kg-'s 

0.02 

0.2 

X 2 

Mobility of 2nd ion 

3.1083 x 10 11 kg~'s 

0.02 

0.2 

£ll 

Self-interaction of 1 st ion 

1.6606 x 10- 5 m 3 mor' 

1 

to- 6 

^22 

Self-interaction of 2nd ion 

1.6606 x 10- 5 m 3 mo/ _1 

1 

to- 6 

^12 

Interaction of the two ions 

1.6606 x 10“ 5 m 3 mo/ _1 

1 

10“ 6 
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If we further remove the ji 3 terms from the ID governing 
equations of the EPNP model, we get the ID equations of the classi¬ 
cal PNP model with the finite size effect (CPNP). In the following, we 
will compare these three models explicitly in ID. First, we examine 
their linear stability properties. 


small |k|. So, the solution is stable. Analogously, the mode of insta¬ 
bility is absent from the CPNP model, which is a limit of the EPNP 
model, at both long and short waves. 

For intermediate waves, we notice a possible instability if b is 
negative, i.e., 


3.4. Linear stability of the constant state 


a > 0, b < 0. 


(3.14) 


If we assume e 0 = 0 (namely, the system does not have 
a permanent charge present), there exists a constant solution 
of the full model, which is a solution of all limiting models, 
ti\ = n 2 = n°,v = 0,p = 0,4> = 0, where n° is constant such that 

n° = r v 0 - r\n° - r v 2 n° > 0. (3.11) 


In certain parameter regimes, the growth rate a, (given in the 
Appendix) can be negative for a very small |fc|, becomes positive for 
some intermediate values of |fc|, and then turns to negative again at 
large | k|. Assuming 7 , = y 2 = 8 « 1, we obtain the roots of A = 0 
asymptotical. Then, we obtain the cutoff wave numbers. We denote 
c = + Sn + -y- + S&r + ?22 + -y-, then we have 


This inequality is necessary to ensure that the solvent density is greater 
than zero. We perturb this constant solution as follows: 

n, = n° + ee at+ikx n°, n 2 = n° + ee at+ikx n°, 

Vi = £e “ t+ife v?, p = ee at+ikx p°, * = ee “ t+ifet <D°. (3.12) 


A/(A 1 \ 2 (n°) 2 k 2 ) = a + bk 2 + c8k 4 + 8 2 k 6 . 

There are two positive roots of k 2 , corresponding to two cutoff 
wave numbers asymptotically: 

(kf off y = - °^-8 + 0(8 2 ), (k c 2 utoff y = X ° + x, +0(6), (3.15) 


Here e « 1 is a small parameter, a is the growth rate and k is 
the wave number. First, we point out that the velocity components 
v 2 ,v 3 are decoupled from the rest of the system in the linearized 
equations and they do not contribute instability in this problem; so, 
we only consider the coupled system involving the remaining vari¬ 
ables: p°, 4>°, v'j 1 , rf, n%. The linearized eigenvalue problem for the Full 
model is given in the Appendix. The asymptotical analysis in the small 
wave number regime shows that the instability can incur only when 
4, 2 is negative enough. But, 4,, > 0 in the model. So this mode of insta¬ 
bility is absent from the full model and its limits. The system is stable 
for long wave (small wave number) perturbation. It is easy to find 
that the system is also stable for short wave (large wave number) per¬ 
turbation. From the numerical studies, we find that the intermediate 
wave instability appears when 4, 2 is sufficiently large. The analytical 
result of the intermediate wave instability is hard to obtain from the 
full model, but easy from the EPNP model. We thus focus on the linear 
stability of the limiting EPNP model in the following. 

The linearized eigenvalue problem is given in the Appendix. The 
instability condition is A < 0, where 


A = \i\ 2 (n 0 ) 2 jafc 2 + bk 4 + [ 72(^0 + €11 + ^ 

( 1 r v r m W 1 

+ JJ ^ + ^ 2,<8 } < °' ( 3 - 13 ) 

Here, the parameters a, b are defined by 


a = 

N, n° N 2 n° 

b = 

7i+72 . 


a + S 11 + + 24,2 + ^ (rV + r 2 ) (r? + r 2 m ) j +, 

/ r v .rT'\ 1 / rXr 2 \ 1 


+ r^4„ 


r m r v 1 r m r v 

+ gll€22-&~ 1 2 0 2 ' ^2- 

n 3 


Because A,, \ 2 , 7 ,, 7 2 ,4,,. 4 22 are a H positive, so the coefficients of 
k 8 ,k 6 are all positive. It implies that A > 0 for large wave numbers. 
This means that the system is stable for short waves. For long waves 
(small wave numbers), since the parameter a > 0, then A > 0 for 


where xo = c+Vc 2 4b _ — ° We only retain the 

first two terms in the asymptotic roots. The parameter A is neg¬ 
ative when the wave number is between the two cutoff wave 
numbers 0 < k™^ < k < k c 2 moff . The growth rate oq is pos¬ 
itive in this intermediate wave number regime. This instability 
depends strongly on the interaction parameter £j, 2 . the instability 
condition is satisfied for a sufficiently large £j, 2 . We fix the parameter 
N, = N 2 = 1,\, = A 2 = 0.02,7, = 72 = 10- 4 ,4„ = 4 22 = 1, 
n° = 1 in the length scale l 0 = him regime, and vary the param¬ 
eter 4, 2 . When 4, 2 > 2.06, b < 0, the intermediate wave instability 
incurs. In Fig. 1 (a), we show the curves of the two cutoff wave 
numbers as a function of 4, 2 . The smaller cutoff wave number 
is decreasing and the larger cutoff wave number k% Jt ° 2 ' is increas¬ 
ing as 4,2 increases from 2.046. The unstable wave number regime 
(kf off ,kf off ) widens as 4, 2 increases. 

When 8 -> 0 , we have -> J ^2 and k 2 ut0 ^ -> + 00 . The 

unstable wave regime is k > k™ 10 ^. The system is unstable for 
large wave numbers (short waves), which is known as the Hadamard 
instability. Hence, the high order diffusion coefficients 7 ,, 7 2 have 
the effect to suppress the short wave instability. 

This intermediate wave instability is also dependent of the con¬ 
stant state n°. When the interaction parameter 4, 2 = 2.2 is fixed, 
but n° is varying, we find that b is positive for small )i° and nega¬ 
tive for large n°. That means the system is stable for dilute solution 
but unstable for rich solution. We also plot the cutoff wave numbers 
as functions of n° with fixed 4 i 2 = 2.2 in Fig. 1(b). The instabil¬ 
ity appears when n° > 0.87, and the unstable wave number regime 
(k™ tof K k c 2 ut0 ^'j widens as n° increases. 

This intermediate wave instability is a feature of these 
three models. Through a numerical investigation, we confirm 
that this instability property can occur in all three models. 
In the following example (Fig. 2), we use parameter values 
Ni = N 2 = 1, A, = A 2 = 0.02,7, = 72 = 10- 4 .4,, = 4 22 = 1, 
4, 2 = 2.2,n° = lin the length scale lo = 1 nm regime. The 
instability condition (3.14) is satisfied. The two asymptotical cutoff 
wave numbers of EPNP model are = 9.24 and k c 2 t0 ^ = 45.23, 
respectively, when length scale lo = 1 nm. For the three models, 
the relation between the length scale and the growth rate follows 
a simple scaling law: we denote the two length scales as lo'./g 2 *, 
the corresponding growth rates as oc^\oc^\ and the cutoff wave 
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Fig. 1. Cutoff wave numbers as functions of (j 12 with n° 



(b) Fixed \\2 = 2.2 


= 1 and as functions of n° with 4,2 = 2 . 2 . 


numbers as k^\k^ 2 \ respectively. If = K, then the cutoff wave 
number ratio follows ^,! = I( while the growth rate ratio follows 
% , lU = K 2 - 5 . This can be inferred from the definition of time scale 

4 W 

f° = J ~ lg 5 ,(rn 0 ~ lg). The numerical results in Fig. 2 also 
confirm this analysis. 

The analysis and numerical results show that the growth rates 
can be positive in some intermediate wave number regime depicted 
in Fig. 2, instead of near the zero wave number range. In this case, 
the growth rate of the Full model is the smallest while the EPNP 
model’s is the highest. From the linear stability analysis, we notice 
that this instability is associated with a large interaction parameter 
£ 12 , a consequence of the finite size effect. A positive £ 12 means that 
the interaction between different species due to their steric effects is 
repulsive. The analysis and numerical results tell us that the interme¬ 
diate wave instability appears when the repulsive effect is sufficiently 
strong in the three models. This also can be obtained from the interac¬ 
tion free energy density g. The repulsive interaction due to the finite 
size effect is represented by ^ + 2£i 2 nin 2 + § 22 ^ 2 )' which 

\ 2 

c/fnn 1 + -%=n 2 J + 

When £ 12 is sufficiently large, £11 £22 - £f 2 < 0. this quadratic form 
is hyperbolic type without lower bound. In the next nonlinear sim¬ 
ulations, we only consider the cases § 11^22 - £12 > 0. with out the 
intermediate wave instability. 



can be rewritten as 


M 

2 


( 


3.4.1. Discussion on the finite size effect 

The hard sphere repulsion characterizes the finite-size effect of 
ions, witch keeps ions apart. The free energy density due to the finite- 
size effect is 


J K'(x-y)G({n i }^ 1 (x),{n i }SL 1 (y))dy = J X X - yp ni ^ n J^ dy - 


i=ij=i 


(3.16) 


where a t and a, are the radii of ion i and j, and £y is the energy cou¬ 
pling constant between ion i and j. Thus, in the free energy function, 
we have the convolution integral with the following form 

/ / j^fl2 ni(x)n J (y)dydx - (3 ’ 17) 


We can approximate the above convolution integral by truncating 
the kernel _^ 12 with the cutoff length 6. As discussed in the paper 
[66], when the cutoff length 5 goes to zero, this convolution integral 
can be approximated by the integral 


S s 


n iW n j(y)dx, 


(3.18) 




Fig. 2. The growth rates of the full model and the two PNP models with length scale / 0 = 1,10, lOOnm, respectively, in the parameter regime of intermediate wave number 
instability. The values of growth rates increase as the length scale l 0 increases. The two cutoff wave numbers of the EPNP model in the length scale l 0 = 1 nm regime are 9.24 and 
45.23, respectively. The full model is more stable than the other two models in this regime. 
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with Ss « 6 12+d , where d is the dimension. The free energy density 
due to the finite-size effect can be written as 


y ( Q i + a j ) u S ti n i (x)nj(x) = k B T ^. =1 yn.-n,-, (3.19) 

with + dj) 12 S fi . We add the conformational entropy 

in terms of the derivative form to compensate for the approxi¬ 
mation error, then the energy density for the finite-size effect is 
approximated by 


g = i< b t 


V 1 ' 


ii[ 

-'ij= 1 2 


mi + X i= , j 'I Vn 'ii 


(3.20) 


where y,- is a small parameter, witch can be zero. In the paper [66], 
the following e,j values for the cross hard-sphere potential terms for 
some familiar ions (Na + ,C/“,Ca 2+ ) are used: 


£wa,Na : £a,ci ■ £ca,ca ■ £ Na.ci ■ £ Na.ca '■ £c/,ca — 1:1:1: 0.955 : 1 : 0.961. 

(3.21) 

Also in the paper [66], the ratios of the interaction coefficients are 
given for some familiar ions (Na + , Cl ~, Ca 2+ ) as follows 

%Na,Na ■ &I.CI : tca.Ca ■ £,Na,CI ■ § Na.Ca ■ &l,Ca 

= 1 : 2280 : 1.64 : 42.2 : 0.642 : 50.4. (3.22) 

It is easy to verify that ^ 11^22 - £12 >0 for two of the three ions. 
For the familiar ions, the interaction coefficients ^ are in the stable 
regime. That is the reason we only consider the stable cases in the 
nonlinear simulations next. 

3.5. Nonlinear dynamics 

We next explore nonlinear dynamics of the models in the linearly 
stable regime. We use the characteristic length scale l 0 = him 
and set the domain as x e [0,10]. The values of the interaction 
parameters are chosen as fj n = £j 22 = 1,§ 12 = 0.8, satisfying 
£11 £22 - > 0 - We a * so set diffusion coefficients Ti = y 2 = 0- 

The boundary conditions for the number densities ni,n 2 are no-flux 
boundary conditions (2.25); for the velocity iq, the boundary con¬ 
ditions are set at Vi| x =o,io = 0, and for the electric potential 4>, 
they are set at <t >| x=0 = 0,<t>| x= i 0 = <t>o, where <J> 0 is the elec¬ 
tric potential at the right boundary x = 10. We set <h 0 = 1 
in the following simulations. The initial conditions are given by 


n\ = n 2 = n° = l.Vj = 0 ,p = 0,<t> = <f> 0 x/10. The given exter¬ 
nal electric potential is cf> e = rb 0 x/10. The dimensionless mobilities 
are given as A.i = \ 2 = 0.02. We compute the ionic number den¬ 
sities using the Full model, the EPNP model and the CPNP model, 
respectively. 

Fig. 3 depicts the final steady states of the Full model 
and the EPNP model, and the difference between them, where 
§11 = £22 = l'£i 2 = 0-8 in the stable regime. The states 

of the number densities are almost identical in the middle of the 
domain, while the visible differences appear near the two bound¬ 
aries. Because the electric potential is positive at the right boundary 
and zero at the left boundary, some negative ions gather at the right 
side while positive ions gather at the left side due to the Coulomb 
force, forming two visible boundary layers. 

As shown in Fig. 3, the density differences between the two mod¬ 
els are about 0(1) x 10~ 2 near the boundaries. As a conclusion, 
the compressibility of the flow in the full model plays relatively 
important role, it impacts the aggregation effect of the ions near the 
boundaries. 

In the above example, the density ratio is chosen as p 3 : p 1 : 
p 2 = 1 : 0.5 : 2. By halving density pj and doubling density 
p 2 to increase the density differences, we reset the density ratio 
as p 3 : pj : p 2 = 1 : 0.25 : 4 and p 3 : p, : p 2 = 1 : 
0.125 : 8, while maintaining the volume ratio unchanged at v 3 : 
Vi : v 2 = 1 : 2 : 1, then the dimensionless parameters are 
Ri = 0.0375, = -0.075 and R, = 0.04375, R 2 = -0.175, 

respectively. In these cases, the size differences of the three compo¬ 
nents become larger. As shown in Fig. 4, the differences between the 
Full model and the EPNP model become larger as the size differences 
become larger. When we halve the density pj and double the density 
p 2 , the absolute maximum difference between the Full model and the 
EPNP model is almost doubled. As a result, when the size differences 
between the components are enlarged, the parameters Rj,R 2 are no 
longer small so that the compressibility of the flow can no longer be 
neglected. 

We also compare the EPNP and the CPNP model in the stable 
regime with = £22 = l-£i 2 = 0.8 in Fig. 5. The density ratio 
P3 : Pi : P2 = 1 : 0.5 : 2 is used. The values of the param¬ 
eters are set at r™ = 1 ,r™ = 2. The number density differences 
between the two models are about 0(1) x 10 -2 near the bound¬ 
aries. The differences near the boundaries in n 2 are bigger than that 
in Hi. The reason is that in the CPNP model, the term r["p 3 is dropped 
in the n ; transport equation. In this example, r^ 1 > rj", so the dif¬ 
ferences near the boundaries of n 2 are bigger. Consequently, the 
solvent’s chemical potential p 3 in the EPNP model plays a relatively 
important role, it impacts the aggregation effect of the ions near the 
boundaries. 





(a) Number density 


(b) Differences between the (c) Electric potential <J> 
models 


(d) Energy F(t) - F(0) 


Fig. 3. Steady states of the ionic densities and the electric potential of the Full and the EPNP models with = £ 2 2 = l.di 2 = 0.8 in the stable regime, respectively. The 
differences appear near the boundary with absolute maximum difference 0.011. The curve of energy difference F(t) - F[ 0) is plotted with respect to time. The total free energy 
F(t) decays to a constant when the final steady state is obtained. 
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(a) p 3 : pi : p2 = 1 : 0.5 : 2 


(b) P3 : Pi : p2 = 1 : 0.25 :4 


(c) P3 : Pi : P 2 = 1 : 0.125 : 8 


Fig. 4. Differences between the Full and the EPNP model with = £22 = l,£i 2 = 0.8 in the stable regime, (a) The density ratio is p 3 : Pi : p 2 = 1 : 0.5 : 2 and 
Ri = 0.025,R 2 = -0.025, the absolute maximum difference is about 0.011; (b) the density ratio is p 3 : p, : p 2 = 1 : 0.25 : 4 and R, = 0.0375,R 2 = -0.075, the absolute 
maximum difference is about 0.022; (c) the density ratio is p 3 : Pi : p 2 = 1 : 0.125 : 8 and Ri = 0.04375,R 2 = -0.175, the absolute maximum difference is about 0.044. 
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(a) Number density (b) Difference between the models 

Fig. 5. Steady states of the ionic densities and electric potential of the EPNP and the CPNP model with = | 22 = 1,| 12 
the boundary layers near both ends of the domain. 


x 

(c) Electric potential <l> 

0.8 in the stable regime. The difference appears at 


Next, we consider ionic concentrations without the finite size 
effect and compare them with ionic concentrations with the finite 
size effect using the classical PNP model. In the following, the 
OPNP means the classical PNP model without finite size effects (i.e., 
£11 = £22 = £12 = 0 and 7i = 72 = 0.) The differences also 
appear in the areas near the two boundaries. As shown in Fig. 6, the 
differences of the ionic density can reach up to 0(1) x 10 -1 . The 
finite size effect plays an important role in the system, it impacts the 
aggregation effect of the ions near the boundaries, as studied in the 
papers [26, 27, 33, 66], 


Based on our numerical investigations and the linear analysis, we 
conclude that the ID steady states of the number densities are nearly 
identical in the middle of the domain in all three models in the stable 
regime. The differences lie in the areas near the boundaries. The com¬ 
pressibility of the flow, the chemical potential of the solvent and the 
finite size effect are three main reasons that lead to the differences. 
So, our quasi-incompressible model (the full model) seems to be 
more reasonable because the mass and momentum conservation 
laws are preserved in the model while the other models don’t respect 
the two fundamental physical conservation laws. 



x 

(a) Number density 


x 

(b) Difference between the models 


X 

(c) Electric potential O 


Fig. 6. Steady states of ionic densities and the electric potential of the CPNP and the OPNP model. In the CPNP model = £ 2 2 = E£i 2 = 0.8 in the stable regime. The 
difference appears at the boundary layers near both ends of the domain. 
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Further investigations in higher dimensions is necessary to eval¬ 
uate the difference among the models, which will be conducted in a 
sequel. 


4. Conclusion 

We have developed systematically a set of quasi-incompressible 
theories for ionic fluids of multiple species that respect not only 
momentum conservation but also mass conservation at the presence 
of the ionic species. The previous PNP type models are approxi¬ 
mations of the more general theories when densities of different 
ionic species are distinct. In these theories, we consider the entropic 
contribution from each ionic species together with the ion-ion inter¬ 
action due to the finite size effect. The limiting cases include the 
extended PNP, the classical PNP with the finite size effect, and the 
classical PNP model without the finite size effect. At the length scale 
larger than hundreds of nanometers, all models agree with the clas¬ 
sical PNP model very well. At the length scale in a few nanometers, 
the models can predict quite different stability behavior for homo¬ 
geneous equilibrium states. In nonlinear dynamics, the ionic number 
densities are nearly identical in the middle of the domain, but the 
differences lie in the areas near the boundaries. Apparently, three 
main factors in the compressibility of the flow, the chemical poten¬ 
tial of the solvent and the finite size effect of the ions can lead to the 
discrepancy in model predictions. We tend to believe that the new 
model is more accurate since it obeys the two fundamental physical 
conservation laws in mass and linear momentum while the others 
don’t. 


Acknowledgment 

Xiaogang Yang’s work is supported by the Scientific Research 
Fund of Wuhan Institute of Technology through grant K201741; 
Jun Li’s work is partially supported by NSF of China through 
a grant (NSFC-11301287); Qi Wang is partially supported by 
NSF through awards DMS-1200487 and DMS-1517347 as well 
as grants from NSFC # 11571032, # 91630207 and NSAF- 
U1530401. 


Appendix A 

The linearized eigenvalue problem for the Full model is formu¬ 
lated as follows, 


//°° 

0 0 
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0 0 p k a 


Ricr 
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0 0 0 a 

0 0 0 0 
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R 2 \ 2 n°k 2 -A 2 n°k 2 
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\ / 


where the parameter values are given by = r|J - rj'n 0 - r v 2 n° and 
p k = 1 - R^n° - R 2 n°. The other components in the matrix are 
defined as follows 

B, = -ikr\ + ikn° (-L + + g 12 ) + ik 3 n°y u 
B 2 = -ikr v 2 + ikn° + & 2 + £ 12 ) + ik 3 n°y 2 , 

Ci = A_i n 0 r™r"^ + A in ° k 2 + Ain°y,k 4 , 

C 2 = A.! n°r^r 2 ^ + Ain°(€ 12 )k 2 , 

Di = A 2 n°r 2 m r^ + A 2 n°($ 12 )k 2 , 

U 3 

D 2 = \ 2 n°r™r %^ + \ 2 n° + € 22 ) l< 2 + \ 2 n°y 2 k 4 . (5.2) 


Although the coefficient matrix is 5 x 5, the characteristic polynomial 
of the coefficient matrix is a third order polynomial of growth rate a, 
which yields three independent eigen-modes. Using an asymptotic 
analysis at small wave numbers |k| <sc 1, the three asymptotic growth 
rates are obtained asymptotically: 


a x = T, k 2 + 0(k 4 ), 

T 2 ±Jq + 


T 3 


£*2,3 = 


T 4 


■ 0(k 2 


(5.3) 


where 


h = 


AiA 2 n° 


nf [ P * n ° (jv^ + + + 522 + 2 ^ 12 ) 


(Ai + A 2 ) n 3 
+ Pk (r\ + t 2 ) (jr™ + r + (R, + R 2 ) n§)], 


T 2 —\i\ 2 Pk(n°) 2 (Ri + R 2 ) + £, 


T 3 — - 4(n°) 2 p k £ (Ai + A 2 ) ^AiR 2 + A 2 R 2 ) , 

T 4 =2 n°p k £ (A,R 2 + A 2 R 2 ) . (5.4) 


Notice that a 23 < 0 for small k due to T 3 < 0. The eigenvalue > 0 
when Ti > 0, 


iv^ + N^0 + ^ 11+fe2+ ^ 12< 


4^r(( r ™ + r 2 m ) + (Ri+R 2 )^)- 

"k Tl 3 

(5.5) 


This is the instability condition for long waves for the Full model. It 
follows from Eq. (3.11) that (r ™ + r 2 ) + (Rj + R 2 )n° > 0. So, the 
instability can incur only when £ n + £ 22 + 2§ 12 is negative enough. 
But, fjy > 0 in the model. So this mode of instability is absent from 
the full model. 

For the EPNP model, only ni,n 2 ,<t> are coupled, the eigenvalue 
problem is given by 



/ - fe2 1 ~l 

| A,n°k 2 C, C 2 
-A 2 n°k 2 D, D 2 



(5.6) 


(5.1) 

















690 


X. Yang, Y. Gong.J. Li et al. / Journal of Molecular Liquids 273 (2019) 677-691 


where 


c, = 


n°rM^ + X i n ° (n^o + *n) fe2 + MV< 4 . 


fc 2 


C 2 = \ 1 n°r7 , rJ- ff +\,n^ 12 k 2 , 

n 3 

Di = \ 2 n°ifr^+\ 2 n 0 $ 12 k 2 , 

U 3 

= \-,n°r™rxt + \,n° (—L 


D 2 = A 2 n°r?r^ + \ 2 n° + § 22 ) fc 2 + A 2 n°y 2 k 4 . (5.7) 


Eliminating <T>°, the system reduces to 


a0\ /C,+A 1 n°| C 2 -A in °I 

0 a 1 + I D, - A 2 n°i D 2 + A 2 n°| 


n° | = 0. 


The characteristic polynomial is quadratic and given by 

il' 


Ci + D 2 + (Ai + A 2 ) n° - 


a + A = 0, 


,1 


(5.8) 


A = C iD 2 -C 2 D, +[A 2 (Ci +C 2 ) + A, (D, + D 2 )] n 0 -. (5.9) 


The two growth rates are given by 

(c, + D 2 + (\, + A 2 ) n° I) + ^(c, +D 2 + (Ai + A 2 )n°|) 


a\ = - 2A 


-4A 


(c, + D 2 + (\, + A 2 ) n° I) + J (c, + D 2 + + \ 2 )n°±y - 4A 

(5.10) 


Notice that C\ > 0,D 2 > 0. So, Re(a 2 ) < 0 and Re(a i) can be positive 
only if A < 0. 
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